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Abstract We derive an expression for the mutual gravitational force and 
torque of two bodies having arbitrary shapes and mass distributions, as an 
expansion in power series of their products of inertia and of the relative 



coordinates of their centres of mass. The absolute convergence of all the power 
series developed is rigorously demonstrated. The absence of transcendental 
functions makes this formalism suitable for fast numerical applications. The 
^ i ■ products of inertia used here are directly related to the spherical harmonics 

coefficients, and we provide a detailed analysis of this relationship. 
■ 
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1 Introduction 



The mutual gravitational potential of two physical bodies contains the infor- 
mation necessary and sufficient to derive their mutual force and torque, and 
thus completely describe their dynamical evolution. However, in general the 
direct derivation of the mutual force and torque from the mutual potential 
can pose substantial mathematical difficulties. In the literature, these dif- 
ficulties have often been alleviated by making assumptions on the physical 
properties of the bodies or on their dynamical configuration, but of course 
this comes at the cost of the loss of generality. 

In p] the mutual gravitational potential of N solid bodies is expressed 
using the spherical harmonics coefficients of each body, but in a formalism 
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based on complex variables for the sake of a compact notation. The expres- 
sions make use of transcendental functions, that are computationally expen- 
sive and thus not prone to fast numerical applications. No explicit expression 
for the mutual force and torque is provided. 

A different approach is followed in [5] and [2], where the shape of each 
body is described by a polyhedron, and the spatial resolution of the mass 
density distribution is limited by the number of simplices. Transcendental 
functions are avoided, at the benefit of fast numerical applications, and the 
limitations imposed by the choice to model the shape of each body with 
a polyhedron can be made arbitrarily small by increasing the number of 
vertexes, at the expense of computational efficiency 

A basic but powerful approach is used in 0], where the mutual gravita- 
tional potential is expressed in a formalism based on cartesian coordinates 
and products of inertia, with no assumptions on the shape and mass density 
of the bodies, and avoiding transcendental functions. This is the most generic 
formalism we could find in the literature. 

In this work, we choose an approach that allows to treat the figure-figure 
interaction of rigid bodies in its most generic form (arbitrary shape and mass 
distribution). Given our previous analysis, the natural choice is to follow [4] 
and determine explicit expressions for the mutual gravitational force and 
torque. We also improve on [4] by demonstrating the absolute convergence 
of all the power series developed, and by finding an explicit expression for 
timnLMN, the main coefficient of the formalism, that was only defined recur- 
sively in [4]. 

The plan of the manuscript is the following: in Sec. j2|) we introduce the 
main formalism that will be used throughout this work; in Sec. ((3]) we derive 
the complete expression for the mutual gravitational force and torque; in 
Sec. ((4]) we demonstrate the absolute convergence of all the series developed. 
Two appendix sections further enrich this manuscript: in Sec. (|A"| we de- 
rive the conversion formulas between the products of inertia coefficients and 
the spherical harmonics coefficients; in Sec. jB]) we provide formulas for the 
products of inertia relative to bodies with a regular shape and an uniform 
density. 



2 Mutual Gravitational Potential 

In this section we introduce the formalism developed in [4 , with some modifi- 
cations and extensions. Consider two bodies 6=1,2 with centers of mass Ob 
and volumes Vb- With respect to the inertial reference system, Oi = [0, 0, 0] T 
is chosen to be at the origin, while O2 is located at [£,77, C] T , and the po- 
sition of a generic point Pb within Vb is Pi = [Xi,Yx, Zi] T and P2 = 
[£ + ^2, fj + Y2, C+ Z 2 ] T . The symbol " T " denotes the transpose operator. In 
the body-fixed reference systems, we have Oi = [0,0, 0] T , Pi = [x\, yi, Zi] r , 
and O2 = [0,0, 0] T , P2 = [x 2 , J/2, 22] 1 ". A vector in the inertial frame can 
be transformed into corresponding vector in the body-fixed frame through a 
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rotation Qi, 



hb fob hb 
mib m 2 b m 3b 
nu n 2 b n 3 b 



such that [xb, y b , Zb] T = Qb[Xb, Yb, Zb] T - Similarly, we define the vector R = 
2 Oi = C] T in the inertial frame, that transforms as {£b,Vb,(b] T = 
Qb[L V, (V, so that R 2 = £ 2 + n 2 + Q 2 = $ + f]\ + Q. The inverse of the 
distance p between Pi and P2 can now be expressed using Eq. (6") in [3]: 
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(1) 



where 



«5 = «i + U + H + h 

35 = h + h - U + h ~ j a - 

h = fcl + «2 - «3 + J2 ~ J3 



fc 3 

■k 2 



■*3 



and i 5 + j 5 + fc 5 = i x + i 2 + ji - 
defined by the recursive relation 



h +k\ + k 2 . The coefficient ti mnLM N is 



timnLMN = { {21 — Cl m n)t(l-l)mn(L-l)MN + 

+ (I — — Ci mn )t(i_ 2 )mnLMN} ($10 ~ 1) + 

+ {(2m — Cjmn)*J(m-l)ni(M-l)JV + 

+ (m - l)(m - Ci„ m )i;( m _2)nLAfAf} (5m0 - 1) 

+ {{2 n — Clmn)tl m (n-l)LM{N-l) + 

+ (n— l)(n - Cl mn )tlrn(n-2)LM n} (#n0 ~ 1) 



(2) 



with Cimn = 1/ (3 — Sio — S m o — 5 n o)i an d % the Kronecker delta function. 
The initial value is ioooooo = lj and the coefficient is zero when any one 
of (I + L), (m + M), or (n + N) is an odd integer, or when any of the six 
indices assumes a negative value. A table containing the coefficients up to 
I + ni + n < 7 is available in [3] , and in Eq. (|13[) we derive an explicit formula 
for timnLMN- Finally, the coefficients l x , m x , ■ ■ ■ , n z arc defined by: 



lx fo 

m x m y m z 

Tlrn Tin Tlv 



QiQl 



(3) 
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where the matrix Qi clearly represents the relative rotation between the 
two bodies. 

The mutual gravitational potential U±2 can now be written as 

U 12 =G [ [ M*l,l/l,*l)fr(* a ,W,*2) d y id y 2 (4) 

JVi Jv 2 P 

where G is the gravitational constant and 5t,(x b , y b , z b ) is the mass density of 
the body b. By substituting Eq. |T]) in Eq. (j4|) , we obtain integrals defined as 
generalized products of inertia in 0] , of the form 

J J J? 



Mb,ijk — / x l b y 3 b z b Sb(xb,yb,z b ) dVb (5) 



and since these integrals involve only coordinates relative to the body b in 
the body-fixed reference frame, they play the same role of the {Ci m ,Si m } 
spherical harmonics coefficients (see Sec. (©), and need to be computed 
only once for a body with constant shape and mass distribution. In this 
manuscript we prefer to deal with adimensional coefficients, so we introduce 
the normalized products of inertia j\f b ,ijk defined by 



l+J+k S b (x b ,y b ,z b ) dV b 

Jvb,ijk — ? — 1T7Tk~n \° 

/ 5 b { XbM )W b r^ k M b , 000 
Jv b 

where ro is an arbitrary normalization radius, and A^h.ooo = M. b is simply 
the mass of the body b. The final expression for the mutual gravitational 
potential is then: 

A A KA oo oo oo oo oo oo 

^ = G M 1 M 1 E E E E E E x 

ii=0ji=0 fci=0 i 2 =0 j 2 =0 k 2 =0 

M.iu'ifei J^2,i 2 j 2 k 2 fr yi+n+ki+i2+j2+k2 

X 1 



h. 32 ka i 3 h k 3 / ■ . , . . , , \ /•->/•>/ , 
^ >p y> y> ( '->\ ( ji\ / '-'2 \ / \ / .)■■; \ ( '>';: 

13=0 J3=0 fe 3 =0 i 4 =0 j 4 = fe 4 =0 



v*3/ \h) \kz) \ij \jj 



x ll i ml 3 - u nl 2 - l ni i ml 3 - H nl 2 - 3 H k z i m k /- ki n k z 2 - k3 x 

*s JB k 5 £L^MqN 

><EEE (7) 

_L=0 M=0 JV=0 

The infinite series appearing in these expressions are typically truncated in 
real- world applications. As we show in Sec. (|A"|). the order of the A/yfe coeffi- 
cients is given by i + j + fc. The level of accuracy required by the particular 
problem being investigated drives the order of the coefficients to be used. 
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3 Mutual Gravitational Force and Torque 

The expression 

-G8 1 {P 1 )8 2 {P 2 )V{\/p) dVidV 2 

represents the force acting on Pi due to the gravitational interaction between 
two infinitesimal masses with volume dVf, and densities Sb(Pb) located at p,, 
with positive direction P 2 — P\. In order to obtain the total force acting on 
the barycenter of body 1, we need an expression for the gradient V(l/p) with 
components: 



V(l/p) 



d d d 
<9£i ' dr]! ' d(x 



1 T 



The only part of l/p affected by the gradient is tfrj^Ci /R p+1 (see Eq. (D), 
where we have defined p = ii + ji + ki+i 2 + j 2 + k 2 + L + M + N for brevity, 
and where R 2 = £ 2 + rf + ( 2 = Ci + ^i + Ci - Explicitly the effect of the 
gradient on the affected part of l/p is: 



d/di]\ 

d/dCi 



iv/Ci 



i? 2 



Ci 



The total force F\ 2 acting on body 1 is determined by integrating over the 
volume of the two bodies, obtaining: 



oc oc oc oo oo oo 



»1 =0.31=0 fei=0 i 2 =0j 2 =0 k 2 =a 



R 
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(8) 



In a similar fashion, the torque acting on the barycenter of body 1 due to 
the mutual force between the same two infinitesimal masses is given by the 
vector product 



-G,5 1 (P 1 )^(P 2 )[x 1 ,y 1 ,z 1 ] T AV(l/p) dlW 2 
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and the total torque t\% acting on body 1 is therefore: 
,MxM.2 



oo oc oo oo 



ri2 = q'_i±i_i± e E E E E E (-i) ii+ - +fei 

ii=0 ji=0 fc 1= 0i 2 =0 j 2 =0 fe 2 =0 



R 

1 Mi,i a j a ka 



ii!ji!fci! i^-h^v- 



To 



! Vi? 



il+jl+ki+i 2 +j2+k 2 



xEEEEEE 

»3=0 J 3 =0 fc 3 =0 14=0 J4=0 fc 4 =0 



J 5 



At, 



EEE' 



L=0 M=0 JV=0 

•A/i,ii, 3 'i,fei+l^l - M,ii + lji,fciCl 
_M,i 1 +l,ji,fe 1 ??l - ■A/l,ti,j 1 +l,fcifl_ 

M,»iA+i,fci-^/Ci - ^i,ii,jiM+i M /vi 

M.ii ji,fci+li/^l - ■A/l,ii+l,ji,fei-^/Cl 



i? 2 



(9) 



Both Eq. (|8I9[) are expressed in the reference system relative to body 1. It is 
important to stress that the only singularity in Eq. (|8l9p is still in i? — > 0, as 
a possibly suspicious term of the type will turn out to be regular when 
each equation is expanded and the term becomes of the type L£f -1 . 



4 Formal Convergence 

In this section, we analyze the convergence properties of the power series 
expansions contained in this manuscript. We start with a basic relation de- 
veloped in section 202 of [3J, for the vectors p = R r: 

i 1 00 

p » (ooflA) (10) 

F 71=0 



with r = ^/ x 1 + y 2 + z 2 , P„(cos A) the Legendre polynomial of cos A, and A 
the angle between the two vectors R and r. The power series in Eq. (flQ|) 
converges absolutely, because |P n (cosA)| < 1 for any integer n > and for 
any angle A, and: 

£j©"*<-H*£(i)"-( i -5r 

n=0 n=0 

for r < R, where both r and R are positive. This expression for l/p can be 
further expanded while conserving the same convergence characteristics, and 
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section 204 of |3] gives: 

oo L"/2J / s T i„A 7 h 

r E£ - <2 "- 2s - 1)!! E ^* 

X ^ a!/3! 7 ! (i - 2a)!(j - 2/3)!(fc - 2 7 )! i?"- 2s 

If we substitute r = [x, y, z] T with T2 r± = [xz — x%, yi — y%, zi — zi] T and 
then apply the binomial theorem, we obtain: 

oo L«/2J , 1 x s 

!4EE H) (2n-2 S -l)!!x 

n=0 s=0 



,a Kb J \cj R n 

i,j,k a=0 b=0 c=0 \ / \ / V / 
i+j-j-fc— n 

x V _J 1 ^2^-2^-2, 

a!/3! 7 ! (<-2a)!(j-2j9)!(fc-27)! i?"" 2s 1 J 

ct+/3+7=s 

In order to compare Eq. (jTTJ) with Eq. |T]), we assume for a moment no 
relative rotation between the two body- fixed reference frames in Eq. (p}, so 
that the product QiQj is equal to the identity matrix. This leads to the 
following simplified expression for Eq. fT]): 



1 , 00 00 00 00 00 00 „il,A r fel„l2,j2„fc2 

i = _!_ ,il P\^\^\^\^\^ ( lV'l + 71+fcl X l gl Z l X 2 f/2 Z 2 

ii=0 ji=0 fei=0 »2=0 j 2 =0 fc 2 =0 



*5 is A15 ^L^M^N 
X E E E ^SJB k$LMN Rh +n +kl + i2 l h + k 2 +L+M+N ( 12 ) 

l=o M=o n=o 

where we have now i§ = i± + 12, 35 = ji + J2, &5 = k\ + k%. The expressions 
in Eq. (fTTj) and Eq. (fT2f can now be compared: the fact that both equations 
represent the same quantity 1/p is not guarantee that they have the same 
convergence characteristics. But if we can demonstrate that one equation 
is nothing more than an algebraic manipulation of the other, then the same 
convergence of both equations is guaranteed. This is the case for Eq. (|11[) and 
Eq. (|12j) , because we can obtain an explicit expression for the coefficient that 
was originally defined recursively in Eq. @ by comparing the two equations: 

, , un -2s ( lY ilj'IAr! (2n-2«-l 



2 J aipij] (i - 2a)l(j - 2f3)\(k - 2j)\ 

(13) 



(_ 1) (^ t +5 LAJW )/2 i]jW ( 3ijk + S LMN - 

2(sijk-s L MN)/2 a \p\j\ L\M\N\ 
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where n — Syfc = i+j + k, Slmn = L + M + N, 2s = n — Slmn, L — i — 2a, 
M = j — 2(3, N = k — 2j, s = a+f3+-f. This remarkable result, expressing the 
key coefficient of Eq. (|12[) using coefficients in Eq. (fTTj) . clearly demonstrates 
that the two series are identical, and their different appearance is only due 
to elementary algebraic manipulation. We conclude that Eq. (fT2|) converges 
absolutely for |r2 — < R. In the general case, when a relative rotation 
between the two body-fixed reference frames is present, the expression for 
1/p in Eq. (Q]) still converges, as we didn't assume any particular value for 
the variables Xb, yb, Zb and £, 77, £. But in this case the convergence condi- 
tion becomes \r2 + ri| < R, the worst case for the non-rotated convergence 
condition. 

The expression for the gravitational potential [7 12 in Eq. ([7]), as defined 
in Eq. (j4|, converges where the expression for 1/p converges, that is, in every 
point outside the two bounding spheres relative to the two bodies, as long as 
the two spheres have no points in common. 

Finally, the expressions for the force F12 in Eq. ([5]) and the torque T12 
in Eq. (|9|) both involve the gradient of 1/p, that can now be derived from 
Eq. (fTOf using the relation (1 — x 2 )Pl l (x) = nP n -i(x) — nxP n (x) and the 
expansion cos A = (x£ + yr/ + zQ)/(rR). The resulting expression for the 
gradient of 1/p is composed by three infinite series that can be generalized 
by: 



00 

/TV 



J2{r) (n + p)Pn- q (cosX) 

with p and q integers equal to or 1. This series converges absolutely: 

(l-p)rR+pR 2 



S I (l) " (n + P)Pn - q (C ° S A) I " S (^) " (n + P) 



n=0 



(R-rf 



for r < R, where both r and R are positive. This completes the demonstration 
that both Eq. © and Eq. ([5]) converge absolutely. 



5 Conclusions 

This manuscript provides the equations necessary to apply the formalism ini- 
tially developed by [4] to dynamical problems involving the gravitational in- 
teraction of an arbitrary number of bodies with arbitrary shape and mass dis- 
tribution. The convergence requirements are that any two bounding spheres 
relative to two bodies never intersect, and this is analogous to formalisms 
based on a spherical harmonics expansion of the gravitational potential. The 
dynamical evolution of an arbitrary number of solid bodies with arbitrary 
shape and mass distribution is completely determined by the expressions for 
the mutual force and torque provided in Sec. ((3|). The products of inertia 
relative to each body can be obtained via direct integration (Eq. (|5l6p ), and 
in special cases, for homogeneous bodies with the shape of a box, a cylin- 
der, or a triaxial ellipsoid, the products of inertia have been computed up to 
any order and the results are presented in Sec. (jBJ. The connection between 
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spherical harmonics coefficients and products of inertia coefficients is ana- 
lyzed in great detail in Sec. ([X]). that provides explicit relations between the 
two. Finally, in Sec. (|4]) we demonstrate the absolute convergence of all the 
power series developed in the manuscript, and in Eq. (|13p derive an explicit 
expression for the ti mn LMN coefficient, that plays a fundamental role in this 
formalism. 



A Conversion Between j\fijk and {Ci m , Si m } Coefficients 

In this section, we derive the conversion formulas between the spherical harmon- 
ics coefficients {C; m ,S; m } and Nijk, for a body having arbitrary shape and mass 
distribution. The spherical harmonics coefficients are typically denned by: 



C lm = — / (r/r yPi m (z/r)cos(m<P)5(r) dV (14) 



M 



v 



M 



Sim = -T7 / (r/r yPim(z/r)sm(m<t>)5(r) dV (15) 



v 



where M is the mass of the body, r = \J x 2 + y 2 + z 2 , ro is an arbitrary normal- 
ization radius, P; m (z/r) is the associated Legendre function of z/r, t<m<f) = y/x, 
8(r) is the mass density at r, and V is the body's volume. With the intent to use 
a light notation, the coefficients C; m and Si m defined above do not include any 
normalization factor. Any normalization factor can be easily added at any point in 
this section. 

In order to proceed, we need to transform the expressions in Eq. (|14I15[) to 
make explicit the dependence on the coordinates x, y, z. Consider the expansion 
for the associated Legendre function: 



prw = (-i) m (i-™ 2 r /2 2- i x 

LI/2J 



> > >"'■( ; )\ 21 2 *) ( /-,»~2A + j. • 



l—m — 2k 



where the notation |_aj represents the floor of a, and (a) m is the Pochhammer 
function of a. Since in the literature the associated Legendre function is defined up 
to a factor (-l) m , here we will use the notation Pi m (w) = ( — l) m P l m (w). We use 
expansion of cos(m0) and sin(m0): 

Lm/2J / \ 

cos(m0) = ]T (-l) fc I ™ J sin 2k ((p) cos m ' 2k (<f>) 

k=0 V / 

L(m-1)/2J / \ 

sin(m0) = E (- 1 ) fc f 2k^i ) sm2 ' :+1 (^) cosm_2fc ~ 1 (^) 
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As we have yjl — (z/r) 2 cos(4>) — x/r and -y/1 — (z/r) 2 sm(<j>) = y/r we can now 
write: 



U/2J |m/2J / 

(r/r ) l Pi m (z/r)co S (m4>) = 2~ l Y, E (-1) P+9 

p=0 <j=0 V 



j i7 1 s 



Li/2J L(m-1)/2J / \ / 

(r/r o )'P im (^/r)sin(m0) = 2- ; ^ £ 

p=o 9=0 \ p / V 



n / 21 - 2p\ 

I x 



x ( 2 x I (I - m - 2p + l) m x m -' q - L y' q+L z l - m '^T^/r l 



But before performing the integrals in Eq. fl 141 15 jl we need to expand r 
( x 2 + y 2 + z 2 y as 



and we finally obtain: 



2p 



Li/2J |m/2j /A / 

c, m = 2- i j2 E(-d p+9 J 

P=0 9=0 \ P J \ 

P P—Vg 

x (l-m-2p + l) m J_] 



l\ 2l-2p\ m 



„,. „,. "-'>'» ! (p->'.-^)! 

X A/"m-29 + 2i/ x ,29+2^ H ,i-m — 2^ x -2i^j, (16) 



L!/2J L(m-1)/2J /A / \ / 



X A/ r m _29-l+2^ x ,29+l+2^ H ,!-m-2i/ x -2i/y (17) 



These two equations solve the direct problem of expressing the spherical harmonics 
coefficients {C; m , Si m } as a function of the normalized products of inertia A/ijfc. The 
terms up to the fourth order for C; m (omitting the trivial values Coo = Aooo = 1, 
Cio = Aooi = 0, Cii = 7Vioo = 0, Sn = A/010 = when the center of mass is at 
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the origin) are: 

C20 - 
C21 - 
C22 - 
C30 = 
C31 = 
C32 = 
C33 = 
C40 = 
C41 = 
C42 = 
C43 = 
C44 = 
and for Si m : 



ATooa-1/2 (M020+M200) 
3Moi 

3 (Moo -AT Q 2o) 
Moa -3/2 (Mm + M201) 
6A/i 2 - 3/2 (A/120 + M00) 
15 (A/201 - A/021) 
I5A/300 - 45A/i2o 

A004 - 3 (M02 + M22) + 3/4M 22 o + 3/8 (Mwo + Moo) 
IOA103 - 15/2 (Mm + M01) 
45(A/" 20 2 -M022) + 15/2(Ao4o -A/ioo) 
105A/301 - 315Ai2i 
105(A/4oo +AT 4o) - 630M 20 

•521 = 3Aon 
£22 = 6A110 

5 3 1 = 6M12 - 3/2 (Mao + A/" 2 io) 

532 = 30A/ill 

5 3 3 = 45Mio - 15A/b3o 

5 4 1 = IOA/013 - 15/2 (Mm + Mix) 

5 42 = 90A/xi2 - 15 (Mao +M10) 

5 43 = 315Mn - 105Mo3i 

5 44 = 420 (M10 -Mso) 

In order to solve the inverse problem, that consists of expressing the coefficients 
Mjfc as a function of Ci m and 5; m , we need to solve for Mjfc the linear system: 

Ci m = flmKjk 
Sim = "liMjfc 

is real and two-dimensional (rectan- 



where each of the two matrices 7; 
gular in general): 



and <t^ 



~ Coo " 




Cxo 




_ Clm _ 





,000 ,,100 

7oo 7oo 

,000 ,,100 

7io 7io 



y°°° 100 



ijk 

Too 

ijk 

7io 



7i 





"Moo " 




Moo 







We are allowed to split the solution of the inverse problem into the inversion of two 
separate matrices because each one of the terms Mj'fc is involved in either Eq. (|16[) 
or Eq. (|17[) . A close inspection of these equations in fact reveals that: 

— both Cim and 5; m contain only terms Mjfc where z + j + fe = Z; 

— Cim contains only terms Mjfc where j = 2q + 2^ is an even integer; 

— Sim contains only terms Mjfc where j = 2q + 1 + 2^ is an odd integer. 



If we consider only matrices 7^ and cr^ containing all and only terms of order 
I = i + j + k, we have that 7^ has dimensions (Z + 1) x |_(Z + 2) 2 /4J , while cr^ has 
dimensions Z x |_(Z + 1) 2 /4J. A total of 2Z + 1 equations involving 2Z + 1 spherical 
harmonics coefficients and (Z + 1)(Z + 2)/2 product of inertia coefficients. As it 
is clear, the problem of expressing the Mjfc coefficients in terms of the spherical 
harmonics coefficients is underdetermined, and has infinitely many solutions (if 
there are any). 



ijk 



12 



B Mijk Coefficients for Homogeneous Bodies with Regular Shape 

With the aim to provide more ready-to-use formulas, we now compute the value of 
•A/ijfc for a body having uniform mass distribution and the shape of box, a cylinder, 
and a triaxial ellipsoid. Because of these particular choices, the expression for Nijk 
provided by Eq. ((6} simplifies as follows: 

/ ^*(*,y,*)dV 
M m = ^ = / x v z * dV (18) 



S(x,y,z)dV Vr o 



v 



V 



Furthermore, because we choose the principal axes of the body as reference system, 
and the shape of each body i s a xis-symmetrical in each of the three principal 
directions, the integral in Eq. (|18|l vanishes when any one of i, j, or k is an odd 
integer. So we will only consider expressions for A/yit where i, j, and k are even 
integers. Expression analogous to those presented here are available in [3]. 



B.l Box 

For a cartesian box with dimensions 2a x 26 x 2c respectively along x, y, and z, the 
volume is V = 8abc and we have: 



iij k 

, a be 

JVijk = 



r i+i+k (i + l)(j + l)(fc + l) 



B.2 Cylinder 

For a cylinder with a circular base of radius r parallel to the x-y plane, and height 
2h along the z axis, the volume is V = 2hirr 2 and we have: 



B.3 Triaxial Ellipsoid 

For a triaxial ellipsoid with equation: (x/a) 2 + (y/b) 2 + (z/c) 2 < 1, the volume is 
V — (4/3)irabc and we have: 

i/2 j/2 fe/2 

_ ab J c p_=i 9=1 fr i 

VV ijfc - d +fe ( I+J + fe) / 2+ 2 



'0 



n ( 2u - !) 
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The coefficients generated by Eq. (fl"9|) up to the eighth order are: 



A/000 


= 1 


yvoo2 


= 1/5 (c /r ) 


yvoo4 


/ot: / 4 / 4\ 

= 3/35 (c /r ) 


» r 

A/022 


1 /or /; 2 2 / 4\ 

= 1/35 (b c /r ) 


A/0O6 


= 1/21 (c 6 /^) 


A/o24 


= 1/105 (fe 2 c 4 /r 6 ) 


A222 


= 1/315 (aVcVn!) 


A/0O8 


= 1/33 (c 8 /r a ) 


A/026 


= 1/231 (& 2 c 6 /r 8 ) 


A/()44 


= 1/385 (b A c A /rl) 


A224 


= 1/1155 (a 2 b 2 c 4 /ro) 



Those not listed here can be easily obtained using the symmetry of Eq. (|19[) for 
permutations of ijk. 
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